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We report a theoretical study of the collective optical response of a two-dimensional array of 
nonlinear cavities in the impenetrable photon regime under a strong artificial magnetic field. Taking 
advantage of the non-equilibrium nature of the photon gas, we propose an experimentally viable 
all-optical scheme to generate and detect strongly correlated photon states which are optical analogs 
of the Laughlin states of fractional quantum Hall physics. 



In the latest years, the hydrodynamic properties of 
quantum fluids of light have attracted a strong interest 
from both theoretical and experimental points of view, 
with the demonstration of superfluid flow in degenerate 
quantum gases of dressed photons (the so-called exciton- 
polaritons) in planar semiconductor microcavities [1] and 
the subsequent observation of vortices and solitons in the 
wake of a strong defect [2]. All these experiments were 
performed in a dilute gas regime where a mean-field de- 
scription based on a generalized Gross-Pitaevskii equa- 
tion is accurate [3, 4]. 

Going beyond this regime requires that the underly- 
ing medium show a sufficiently large optical nonlinear- 
ity to induce strong effective interactions between the 
photons. A first step in this direction has been the ob- 
servation of strong antibunching in light emission from 
single- mode cavities in both the visible [5, 6] and mi- 
crowave domains [7] via the so-called photon blockade 
effect. The present experimental challenge is to scale up 
the impenetrable photon regime to arrays of many cou- 
pled cavities, for which theoretical work has anticipated 
the onset of different kinds of strongly correlated pho- 
ton states, from Mott insulators [8], to Tonks- Girardeau 
gases in one- dimensional geometries [9, 10]. 

In the meanwhile, several proposals have appeared for 
generating artificial magnetic fields for neutral quantum 
particles. The Berry phase [11] accumulated by an op- 
tically dressed atom while adiabatically moving in space 
can be described in terms of an artificial gauge field [12]; 
the nucleation of quantized vortices in a dilute Bose- 
Einstein condensate under the effect of such a field was 
observed in [13]. Many authors have speculated on the 
possibility of observing quantum Hall liquid states in 
strongly interacting atomic gases in free space [14] or op- 
tical lattices [15-17]. Very recently, the extension of this 
idea to photons has been theoretically investigated in a 
number of configurations, for instance, arrays of coupled 
optical cavities confining single atoms [18], microwaves 
in circuit-QED devices [19], solid-state photonic devices 
in the infrared or visible spectral range [20-22] . 

The present paper reports a theoretical study of the 
optical response of a coupled cavity system to a coherent 
laser field in a regime where impenetrable photons experi- 
ence a strong artificial magnetic field. Our theoretical de- 



scription is based on a generic model that can be used to 
describe a number of different physical systems, ranging 
from macroscopic optical cavities containing atoms [5] , to 
photonic crystal cavities [6] and circuit QED devices [7] . 
In contrast to prior work [18], we take full advantage 
of the driven-dissipative nature of the photon gas [23] to 
propose an all-optical protocol to generate and character- 
ize strongly correlated photon states which are analogs 
of the fractional quantum Hall (FQH) states of electrons 
in two-dimensional geometries under a strong magnetic 
field [27]. The quantum Hall nature of the generated 
state is assessed in terms of its overlap with the Laughlin 
wave function [28]. We anticipate that detailed informa- 
tion on the microscopic structure of the many-body wave 
function can be experimentally extracted from the field 
quadratures of the secondary emission from the device. 

We consider a two-dimensional square lattice of cou- 
pled cavities under a uniform artificial magnetic field. 
Each cavity is assumed to sustain a single photonic mode, 
to be coupled to its nearest neighbor by photon tunnel- 
ing, and to exhibit a large optical nonlinearity leading to 
strong on-site interactions between photons. The isolated 
system can then be described by the standard single-band 
Bose-Hubbard Hamitonian [10, 22]: 

Ho = Y^ hujoblk -hJ^ Hh^'^'' + h^J2^i(^i - 1)' 

where (hi) is the bosonic creation (annihilation) oper- 
ator for site i and fii = SJS^ is the corresponding number 
operator. The effect of the artificial magnetic field is in- 
cluded via the tunneling phase cpij; J is the tunneling 
strength between nearest neighbor sites (i, j), U is the 
on-site interaction energy, and cjo is the natural cavity 
frequency. 

In the non-interacting case C//J = 0, the spectrum of 
the isolated system Hamiltonian Hq for an infinite lat- 
tice is the Hofstadter butterfly [29] which appears as a 
fractal structure in the energy versus a = (27r)~^ W^^ij 
plane, where the loop sum is performed around a pla- 
quette of the lattice. For a real magnetic field B and a 
lattice spacing a, a would correspond to the number of 
magnetic fiux quanta per plaquette, a = Bo? / {h/ e). In 
the presence of finite on-site repulsive interactions, for 
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filling V = 1/2, the bosonic Laugtilin wave function turns 
out to be the almost exact ground-state wave function in 
the continuum limit a <C 1 [14-16] and to capture the es- 
sential features of the ground state up to a critical value 
of a ^ 0.3 for sufficiently strong interactions U/J>5 
[15]. 

To accurately describe a photonic system, one has to 
account for the finite photon lifetime and for the opti- 
cal pump that is used to continuously replenish the pho- 
ton gas [10]. A coherent pump can be included in the 
model as an additional Hamiltonian term in the form 
HdrWe = Xli + h.c.]. Here we restrict our at- 

tention to the simplest case of a monochromatic pump 
at frequency Up that drives all sites with the same am- 
plitude, Fi{t) = F e~^^p^ . Photon losses at a rate 7 are 
described at the level of the master equation dp/dt = 
— ^[Hq + i^drive, p] + ^[p] with a dissipative term in the 
Lindbladform/:[p] =7^:. [bipbl-{blbip^pblbi) /2] [30]. 
As in previous work [10] the steady-state density matrix 
pss is obtained from a numerical determination of the 
stationary point of the master equation dp/dt = 0. More 
details are given in the Supplemental Material. 

We now present the results of our numerical calcu- 
lations for a 4 X 4 lattice using a basis of states with 
total photon number N = 0,1,2 in the hardcore limit 
U/J = 00, where double occupation of a single site is 
not allowed. Excitation to higher TV states is negligible 
as long as the pump amplitude is well below saturation 
F <C 7. We work in the Landau gauge A = {—By^O) 
with an effective magnetic field strength B such that 
a = 1/4. This choice corresponds to = 4 flux quanta 
through the whole lattice and therefore to a filling frac- 
tion u = N/Na = 1/2 for an TV = 2 state. Periodic 
boundary conditions are assumed [31]. 

The first observable we consider is the total number of 
photons riT = present in the steady state of the 

system, a quantity which in many geometries is propor- 
tional to the total transmitted intensity. This quantity is 
plotted in Fig. 1(a) as a function of the relative pump fre- 
quency Acjp = Up—Uo for three different values of the loss 
rate 7 and a fixed pump amplitude F/7 = 0.1. The main 
feature is a strong peak at Aujp/ J ~ —2.83 corresponding 
to a one-photon transition from the vacuum to the lowest 
one-particle eigenstate of Hq at energy h{uj^^'^ +co'o): the 
position of the resonance is at Aujp = uj'^^\ In addition, 
on the spectrum for the lowest value of 7 we can notice a 
much weaker peak at Aup/J ^ —2.52 corresponding to 
a two-photon transition from vacuum to a two-particle 
eigenstate of Hq of energy h{uj^'^^ -^2ujo)'- in this case, the 
position of the resonance is given by Aujp = ijj^^^ /2 [10]. 

In order to isolate the two-photon peaks from the much 
stronger one-photon background, in Fig. 1(b) we plot- 
ted the ratio P2/P1 where the steady-state probability 
of having = 1,2 particles in the system is defined as 
Pi^2 = Tr[pssni^2] in terms of the projectors Tli^2 onto the 
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FIG. 1: (Color online) (a) Expectation value of the total num- 
ber of particles, (b) Ratio P^/Px of the probability of having 
two particles to one particle in the steady- state, (c) Overlap 
O with the optimized Laughlin wave function. All curves in 
(a-c) are plotted as a function of pump frequency Acjp for 
different values of the loss rate 7/J = 0.002 (blue solid), 0.01 
(red dashed), 0.05 (black dash-dotted). Pump amplitude is 
F/7 = 0.1 for all cases. A 4 x 4 lattice with periodic bound- 
ary conditions is considered in the impenetrable photon limit 
U / J — 00. Vertical solid (dashed) lines indicate the position 
of two- (one-) photon transitions as predicted by the eigen- 
frequencies of Hq. The green thin curves in (b) and (c) refer 
to the non-interacting photon case [// J = for 7/ J = 0.01. 

one- or two-particle subspaces: in this way, clear peaks 
appear at half the frequency of the lowest and the third 
lowest N = 2 eigenstates of Hq . The second lowest N = 2 
eigenstate of Hq does not contribute to this spectrum 
as the matrix element for the corresponding two-photon 
transition appears to be very small for the chosen spatial 
profile of the pump. 

A common procedure in the theory of the quantum 
Hall effect to analyze the physical nature of a state is to 
calculate its overlap with an ansatz wave function, for 
example, a Laughlin wave function. Here we extend this 
approach to a driven-dissipative system: under a weak 
pump condition F <C 7, the component of the density 
matrix on the N = 2 subspace p^^^ = YL2pss^2 involves a 
single pure two-particle state, whose wave function is pro- 
portional to the two-photon amplitude i/jij = Tilpssbibj], 
which then plays the role of a two-photon wave function. 
The quantum Hall nature of the two-photon peaks can 
then be assessed by calculating the overlap of tjjij with the 
pair of generalized Laughlin wave functions ^^^'^^(z^, Zj) 
for our toroidal geometry [32, 33], the explicit form of 
which is given in the Supplemental Material. The opti- 
mized Laughlin wave function is the linear combination 
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of that gives the maximum value O of the overlap 



^=1,2 



(2) 



A plot of O as a function of Acjp is shown in Fig. 1(c). 
For the weakest loss rate 7/J = 0.002, the overlap is 
broadly peaked around the lowest two-photon resonance, 
with a peak value as high as 98.9% (as usual, perfect 
overlap would correspond to O = 100%) and then de- 
cays to zero as the higher two-photon resonance is ap- 
proached: this is a clear signature that the lowest two- 
photon peak indeed corresponds to the transition to a 
two-particle Laughlin state, while the second peak cor- 
responds to some orthogonal excited state. For a more 
realistic loss rate 7/J = 0.05, the maximum overlap is 
somewhat reduced by the reduced selectivity of the co- 
herent drive and the consequent mixing with other states; 
however, its value is still as high as 90.9%. The situation 
is of course completely different in the [// J = case 
of non-interacting photons. In this case, only the one- 
particle peak at Acjp/J ~ —2.83 remains visible in the 
spectrum [thin curve in Fig. 1(b)] and the overlap with 
the Laughlin wave function is dramatically reduced at all 
frequencies [thin curve in Fig. 1(c)]. 
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FIG. 2: (Color online) (a) Sketch of the proposed experimen- 
tal setup to measure the two-photon amplitude V^^j . (b) i = io 
cut of the (normalized) two-photon amplitude '^l^i^j (red) and 
of the optimized Laughlin wave function ^{^Zi^^Zj) (black). 
The two functions are almost indistinguishable. System pa- 
rameters: Acjp/J = -2.8144; 7/J = 0.03; F/7 = 0.1. x,y 
axes refer to the Zj — Xj + iyj coordinate. Reference site io is 
marked by a circle. For each lattice site, the arrow indicates 
the complex amplitude of the wave function. 

In contrast to standard condensed matter systems, 
quantum optical techniques provide experimental access 
to the overlap O with the Laughlin state [34]. In the in- 
frared/visible range, the two-photon amplitude ipij can 
in fact be measured with standard homodyne detection 
techniques [30] by homodyning the secondary emission 
from sites i and j with the coherent pump at Up as 
sketched in Fig. 2. To obtain the different quadratures 
^0? = (^ie"*^^ + 6je^^^)/2 of the field emitted by site z, 
one has simply to adjust the phase delay of the ho- 
modyne beam. As quadrature operators for two different 
sites {i ^ j) commute, they can be simultaneously mea- 
sured and correlated. The full two-photon amplitude i/jij 



is finally reconstructed by combining four separate mea- 
surements: 
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In the microwave domain of circuit-QED systems, the sit- 
uation is even simpler as field quadratures can be directly 
measured using linear amplifiers [35]. A measurement of 
the two-photon amplitude ipij provides an intuitive pic- 
ture of the complexity of two-photon FQH states. A nu- 
merically simulated image of this quantity (actually, a cut 



of ipij for a fixed 



3) is illustrated in Fig. 2(b) and 



compared with the Laughlin wave function; the agree- 
ment is excellent as expected from the high values of the 
overlap mentioned previously. 
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FIG. 3: (Color online) System with hard- wall boundary con- 
ditions and finite interaction strength U / J — 20. (a) Expecta- 
tion value of the total number of particles (black solid curve) 
and P2/P1 ratio (red solid curve) as a function of pump fre- 
quency, (b) Overlap O of the two-photon amplitude ^l^ij with 
the Laughlin wave function. Vertical solid (dashed) lines cor- 
respond to two- (one-) photon transitions, (c) i = io cut of 
the two-photon amplitude ipi^j (red) and of the optimized 
Laughlin wave function (black) for a pump on resonance with 
the two-photon transition at Acjp/J = —2.3120 indicated by 
an arrow in (b); 7/J = 0.01; F/7 = 0.1. 

We conclude the paper by discussing how all the 
previously described phenomena are robust against ex- 
pected experimental imperfections and are not restricted 
to small lattices. Although periodic boundary condi- 
tions may be realized in a photonic system (for instance, 
by suitably connecting sites on opposite sides with opti- 
cal fibers), an experimentally much less demanding op- 
tion is to use hard- wall boundaries. Furthermore, there 
might be structural defects leading to a small inhomo- 
geneity in the lattice potential and interactions between 
the particles might be not strong enough to be in the 
hard-core limit U/J = 00. Transmission spectra for a 
more realistic case of 4 x 4 lattice with finite U/ J = 20 
and an energy offset on the order of 7 on a few lattice 
sites are shown in Fig. 3 together with the corresponding 
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overlap with the optimized Laughhn wave function for 
periodic boundary conditions. In spite of the different 
hard-wall geometry and the disturbing effect of disorder 
and finite interactions, the overlap on resonance with the 
Acjp/J = —2.3120 two-photon transition is as high as 
^ 76%. This fact is intuitively visible in the two-photon 
amplitude plotted in Fig. 3(c). Although there are some 
deviations, the complex structure of the Laughlin wave 
function is still clearly recognizable. 

Calculations for larger systems and more particles can 
be performed with a more sophisticated Monte- Carlo 
wave- function technique [36] which is detailed in the Sup- 
plemental Material and is much more time-consuming. 
For this reason, we did not calculate the whole spectrum, 
but we restricted ourselves to a single value of the pump 
frequency, chosen to be on exact three-photon resonance 
with the lowest three-particle eigenstate of an isolated 
5x5 lattice with a gauge field intensity a = 0.24 so as 
to give again N/N^ = 1/2. For 7/ J = 0.01, F/7 = 0.1, 
averaging over 500 independent Monte Carlo realizations 
yields an overlap of ~ 93% between the three-photon 
amplitude i/jijk = Tilpssbibjbk] and the Laughlin wave 
function, meaning that our all-optical technique is able 
to generate also more complex Laughlin states of more 
than two particles. Of course ijjijk can again be mea- 
sured with the same homodyne techniques by combining 
eight measurements of the form {X^^ X^^ X^^) (cf. the 
Supplemental Material). 

In conclusion, we have shown how FQH states of im- 
penetrable photons can be generated in a small two- 
dimensional lattice in a non-equilibrium setting with sub- 
stantial losses. A particular FQH state is optically se- 
lected by the frequency of the coherent pump and the 
microscopic structure of the wave function can be recon- 
structed from a few homodyne measurements on the sec- 
ondary light emission. We hope our results will trigger 
further work on the intriguing features of FQH physics 
with photonic systems, including fractional statistics and 
topological protection of quantum states. 
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Supplemental Material 



I. NUMERICAL PROCEDURE TO 
DETERMINE THE STEADY-STATE OF THE 
SYSTEM 



A. Super-operator Approach 

In this approach, we first transform the time- 
dependent Hamiltonian Hq + i^drive into a time- 
independent one through a unitary transformation of the 
form U = ex.-p{iujpt^-hi). We then define a super- 
operator S in terms of which the master equation can 
be written as dp/dt = 5p, where p is now a vector con- 
structed from the elements of p. Finding the steady-state 
value pss, which satisfies the equation dpss/dt — Spss = 
0, then reduces to finding the eigenvector of the super- 
operator S with zero eigenvalue. Note that the size of 
the matrix representation of S is x D^^ D being the 
dimension of the physical Hilbert space, which puts a 
severe limit on the system size that could be efficiently 
handled with this method. 



B. Monte-Carlo Wave- Function Approach 

In this method (see Ref. [36] of the main text for 
more details), which can be shown to be equivalent to 
the master equation treatment, the wave function of the 
system is evolved with the non-Hermitian (time- 

independent) Hamiltonian Hq + i^drive — ^^7 ^i/^ and 
at certain random points in time, depending on the occu- 
pation probability of excited states, the system is made 
to spontaneously emit a photon, or in other words to 
make a quantum jump. After running this simulation 
long enough to obtain a sufficient number of statistically 
independent samples of the density matrix \ip{t)){ip{t)\^ 
one can get an estimate of the steady-state density ma- 
trix by averaging over the samples. The different samples 
are obtained from a single quantum jump trajectory with 
sampling intervals on the order of I/7 so that the sam- 
ples can be considered uncorrelated. We checked this 
method against the super-operator method for several 
cases, some of which are presented in the main text, and 
obtained excellent agreement. 



II. LAUGHLIN WAVE FUNCTION ON A 
TORUS 



In an L X L torus geometry, the u = 1/2 Laughlin 
wave function of TV particles in the Landau gauge A = 



{—By^O) has the form 
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where = ^ iyu is the complex coordinate of 
the kih. particle in units of lattice spacing a, Z = 
^^z^, a = Bo?/{h/e) is the number of magnetic ffux 
quanta per plaquette and Ml is the normalization fac- 
tor. The part containing relative coordinates is writ- 
ten in terms of the elliptic ^ functions '^^[d](^k) = 
^^g27rr(n+c) +27r*(n+c) ^ where 71 ruus ovcr ah inte- 
gers. The center-of-mass part is also given by elliptic 
functions: 



1/2 + (iV„ - 2)/4 
-{N^ - 2)/2 



2Z 

T 



2i 



(S2) 



Here, is the number of ffux quanta contained in the 
L X L lattice and the label / = 0, 1 indicates the two 
degenerate ground states at fflling u = N/Na = 1/2. 




FIG. 1: (Color online) Second-order correlation function 
9^'^\'^o,j) for different sets of sites j — {io, ji, . . . ,^5}. Sites 
are ordered with respect to their distance r from the reference 
siteio. (a) {6,10,11,14,15,16}. (b) {6,7,11,8,12,16}. (c) 
{10,6,7,2,3,4}. (d) {10,11,7,12,8,4}. Site enumeration is 
given in the inset of panel (a). Periodic boundary conditions 
and impenetrable photon limit U / J — 00 are assumed. The 
pump is on resonance with the two-photon transition to a 
Laughlin state at /\ujp/J = -2.8144; 7/J = 0.01; F/7 = 0.1. 
Black dashed and red solid lines refer to the Laughlin wave 
functions and to the numerical results, respectively. 
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III. SECOND-ORDER CORRELATION 
FUNCTION 

Numerical predictions for the (here unnormahzed) 
second-order correlation function g^'^\i^j) = (blb^jbjbi) 
are illustrated in Fig. 1 for a 4x4 lattice with a = 1/4 and 
periodic boundary conditions, and a pump on resonance 
with the lowest two-photon peak at Aup/J = —2.8144. 
Specifically, we plot g^'^\io-,j) as a function of the dis- 
tance between site j and a reference site io for different 
sets of sites and we compare it with that obtained from 
the optimized Laughlin wave function which gives the 
highest overlap. The numerical simulation correctly re- 
produces the nontrivial behavior of correlations peculiar 
to FQH states that is visible as peaks and dips in panels 
(c) and (d). 
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FIG. 2: (Color online) i = io cut of the (normalized) two- 
photon amplitude ipi^j for non-interacting photons, U/ J — 0. 
In this case, the cut ij^i^j is proportional to the wave function 
of the lowest one-particle state 0^ at /S.ujp/ J — —2.8284. Sys- 
tem parameters: 7/J = 0.01; F/7 = 0.1. x^y axes refer to 
the Zj = Xj -\-iyj coordinate. Reference site io is marked by a 
circle. For each lattice site, the arrow indicates the complex 
amplitude of the wave function. 



IV. TWO-PHOTON AMPLITUDE FOR 
NON-INTERACTING PHOTONS 

The two-photon amplitude i/jij for non-interacting pho- 
tons with U/J = Ois displayed in Fig. 2. In this case, 
the two-photon amplitude reduces to the tensor square 
tijij = (j)^ 0^ of the wave function (j)^ of the lowest single- 
particle eigenstate with zero momentum along x and does 
not show any special structure in contrast to the Laughlin 
structure shown in the main text for interacting photons. 
For the case of a 4 x 4 lattice with a = 1/4 under investi- 
gation here, each single-particle Hofstadter hand contains 
a single, four-fold degenerate level. The kx = ^ condition 
imposed by the pump profile selects one single state in 
the lowest band and non-resonant excitation of states in 
the upper bands is negligible under the 7/J <C 1 condi- 
tion. 



V. THREE-PHOTON AMPLITUDE 

The A^-photon amplitude V^i^.^i^ = (^ii-.-^i^v) can 
be found by using hi = Xq^ + ^^^^2 obtained from 
the general definition of the field quadrature X^^ = 
(6^6"*"^^ + 6je*'^^)/2. So, it is in principle possible to 
determine this amplitude by doing 2^ separate measure- 
ments. In particular, the three-photon amplitude can be 
written as 



J, — / vi'i') vU) / v(^) vO') v(^) \ 

Wijk - Aq Aq ; - ^Aq A^^2^7r/2/ 

+i{xfx%X^^^) + ). (S3) 



